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We apply an exact diagonalization method to the the infinite-D two- 
band Hubbard model. The method is essentially exact for the calculation of 
thermodynamic properties for all but the smallest frequencies and yields a 
resolution unavailable in Monte Carlo calculations. We establish the insta- 
bility of the normal state with respect to singlet superconductivity at small 
frequencies at small doping, the regime of relevance for High-T c supercon- 
ductors. We also present evidence for the existence of an instability towards 
triplet superconductivity in the large doping regime n ~ 2. 
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The discovery of high-temperature superconductivity has enormously increased the in- 
terest in strongly correlated electron systems. Superconductivity in its various guises has 
been searched for in practically all available models. Of particular interest are models with 
several electronic bands, such as the three-band model of Emery [0 and of Varma, Schmitt- 
Rink and Abrahams 0. These have been proposed as minimal models for the CuO planes in 
the new materials. In 2D, the three-band model has been been actively studied numerically, 
but even the most ambitious Quantum Monte Carlo calculations P] have up to now given 
only inconclusive evidence. 

Given the enormous complexity of the full problem in two or three dimensions it is natural 
to bring to bear the methods of the - now well established - infinite dimensional approach 
on this problem 0]. As for classical statistical mechanics systems, the limit of infinite 
dimensions constitutes a self-consistent mean field theory in which the spatial fluctuations 
are frozen and where the infinite-D many-body system is reduced to a self consistent single 
site problem. It is of considerable interest to clarify whether such models with several bands 
are able to show superconductivity even (or: at least) in large D, i.e. in the absence of 
spatial {e.g. antiferromagnetic) fluctuations. This is the aim of the present paper. 

We investigate in detail the infinite-D two band Hubbard model, which was proposed in 
a recent paper and studied with an extension of the Hirsch-Fye Quantum Monte Carlo 
(QMC) algorithm ||. This method proved to be adequate for the investigation of the normal 
state, and several very interesting metal-to-insulator transitions were identified. However, 
the QMC algorithm did not allow a controlled investigation of the superconducting state or 
of the instability of the normal solution. The main result was that, at the temperatures ac- 
cessible to the QMC method, the superconducting susceptibilities were large in two different 
regimes, the small doping regime (of relevance for high-T c superconductors) and a regime 
with density n close to 2 for intermediate values of the Coulomb interaction. 

In this paper we take up these two regimes [n ~ 1 and n ~ 2). We are interested in 
the normal state exclusively as a starting point for a linear stability analysis and investi- 
gate the regime close to the normal solution. Before presenting our results we will give a 
brief explanation of our exact diagonalization method (cf 0) as applied to the problem of 
superconductivity. Just as for the Hubbard model |7[], this method is far superior to the 
Quantum Monte Carlo method and allows us to decrease by at least an order of magnitude 
the limiting temperature (smallest accessible frequency). 

The (CuO)d-type which is the object of our study is defined by the Hamiltonian || 

H = - J2 tijdtrPja + h.c. + e p PjaPjv + e d d t d iv + U dJ2 n n n< il (1) 

?GDjGP,<t jeP.cr ieD.a i<=D 

where the hopping is scaled as tij ~ l/y/2z (z is the connectivity of the lattice). In eq. ([[]) 
(dcrjPcj) represent two atomic orbitals on different sublattices (D,P) of a bipartite lattice 
with z — > oo which, for the purposes of the present paper, is taken to be the infinitely 
connected Bethe Lattice. The repulsive Coulomb interaction acts only on the ('copper') 
orbital d a , while the ('oxygen') orbital p a is uncorrelated. In the standard Nambu notation, 
= (df, di) (equivalently for ty p ) the d-orbital Green's function can be written as a 2 x 2 
matrix 

D M .-T< W ; H >.(W (2) 
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and the self-consistency equations for the Green's functions are given by || 



D 1 (ico n ) = iu n + {n- e d )a 3 - t 2 pd a 3 P(iuj n )a 3 
V l {iu n ) = iu n + (fj, - e p )cr 3 - t 2 pd a 3 B(iuj n )a 3 



(note that D and D are 2x2 matrices and that D 1 denotes the matrix inverse). 

The Green's function D(iu n ) and the bath Green's function D (iuj n ) are connected by 
the single-site action 



Ssu P = U d f dm d ^r)n dl (r) - f dr f d^^Do" 1 ^ - r')^(r') (4) 
Jo Jo Jo 

This eq. @ is at the base of the Quantum Monte Carlo approach to the solution of the 
d = oo equations, which consists in a decoupling of the interaction term in Ud and a proper 
discretization of the integrals. Alternatively, the Green's functions D(iu n ) and D (iu) n ) may 
be viewed as impurity Green's functions of an effective Anderson model || which, in the 
presence of a superconducting medium, is given by 

n s n s n s 

Ham = e i a ta a iv + Un 11 n n + ^ [VJa^a^ + h.c] + ^[a^a^ + h.c] (5) 

a,l=l a,l=2 l,k=l 

In eq. (^) we have used a compact notation with a site index I, 1 < I < n s , in which 
/ = 1 denotes the impurity. The matrix provides explicit pairing terms between all the 
sites. The matrix is to be taken antisymmetric for the calculation of the triplet sector and 
symmetric in the singlet sector. Of course a strictly self-consistent solution based on the 
correspondence between eq. @ and eq. @ will only be possible in the limit of an infinite 
Anderson model n s = oo. 

The exact diagonalization method consists in approximating the function D (iuj) in eq. 
(0) by the Green's function of a (superconducting) impurity model eq. (|5|) with a finite, and 
even small number n s of sites. Explicitly, Dq^^u)^ 1 = (iuj —H) with H the single particle 
Hamiltonian (U=0) corresponding to eq. (|5|). [In the normal state, in which = 0, the 
Green's function is given by the standard formula G diiu)' 1 = ioj—td—Ud/^—Y^k^ Vk I (^~ 
<k)\- 

In practice, we use a conjugate gradient method to determine the parameters Vk, e^, and 
4>ij minimizing the function 

E=Y, \D (ico n ) - D^ d (ioj n )\/(n max + 1) (6) 

n=0 

Here the uj n denote the smallest Matsubara frequencies. Even at zero temperature we use 
u) n = (2n + l)n/ f3 with (3 a fictitious temperature. This provides a lower frequency cutoff. All 
the results we present do not critically depend on the specific choice of the fitting function, 
and on n max . 

For n s < 6 we are able to diagonalize exactly the Hamiltonian eq. (|^), with a full 
calculation of the spectrum and of the eigenvectors. At temperatures accessible to QMC, 
the Green's functions, susceptibilities, etc can then be compared with the extrapolated 
results of the Monte Carlo algorithm to very high precision. Since we are mainly interested 
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in Green's functions in the zero-temperature limit, we may use the Lanczos method JJ] and 
are then able to increase noticeably the size of the cluster. For the full superconducting 
Hamiltonian eq. (pj) values up to n s ~ 8 may be handled without problems on a work 
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The algorithm which we have just sketched is then iterated to convergence. Once this is 
reached we are able to estimate the quality of the calculation, which is entirely determined 
by the agreement of the two Green's functions Dq and DQ nd . (Note that this is an intrinsic 
criterion). The agreement between the two functions has been found to be amazingly good. 
We insist that the fit of D by DQ nd is the only non-exact part of our procedure and seems 
to be a much less violent approximation than those introduced by the Trotter breakup 
(discretization of the integrals in eq. (f|)) and by the stochastic noise of the Monte Carlo 
procedure. 

An example of the excellent quality of the solution is given by the normal state results 
|TT| at zero temperature presented in fig. 1. Here, the functions Gdo(iw) and G^Q d {iui) are 
displayed. The 'effective' temperature is (3 = 250, Ud = 8, and the density corresponds to 
the lightly doped regime of the two-band model. At n s = 8, the differences are below the 
resolution of the graph with a maximum difference on the order of 10 -2 . At fixed value of 
n max (here taken to be 64), we have noticed a systematic improvement of the quality of fit, 
expressed by a value of E in eq. (|6]) which decreases by a factor of ~ 4 each time one more 
site is added |L2 |. As an illustration of this improvement we show in the insets of fig. 1 



the self-consistent solutions at n s =4,5,6,7, and 8. Above ioo ~ .5 the Green's functions are 
strictly identical. 

We now consider the stability analysis of the normal state solution. A possible way of 
studying this stability is to calculate the pairing susceptibility. An alternative way used 
here is to establish the stability properties of the solution by introducing small terms 0^- in 
the Hamiltonian eq. (H), and following the evolution under subsequent iterations ||. Under 
such conditions, the normal state solutions in fig. 1 very quickly acquire non-zero values of 
F(co), which indicate a superconducting instability. More rigorously, and in order to study 
quantitatively the effects of increasing n s , we may calculate the largest eigenvalue, and the 
corresponding eigenvector of the matrix dF{iu) n+1 / dF{iuj) n close to the normal state, where 
the superscripts on the F's indicate two subsequent iterations of the self-consistency loop. 
This involves a simple rescaling of the <fi at each iteration. We have done such calculations, 
which correspond to the well-known procedure of extracting the largest eigenvalues and 
eigenvectors of a matrix with the 'power method', starting from the parameters in fig. 1, for 
5 < n s < 8. We are able to identify a linear regime at small 0, with the largest eigenvalue 
always of the order \ max ~ 2. The corresponding (rescaled) eigenvectors for n s = 6, 7, 8 
are plotted in fig. 2. Clearly, the agreement between these completely independent curves 
is excellent. We have checked this result in a variety of ways (by changing the effective 
temperature, the precise form of the function used in eq. (0), and the doping). This leads 
us to the conviction that the normal state solution of the d = oo model at small doping is 
indeed unstable with respect to singlet superconductivity. 

We have also studied the point investigated previously ||, i.e. values of the physical 
parameters corresponding to a total density of n ~ 2, where the Hubbard interaction is 
just large enough to create a large overlap between the upper Hubbard band of the <i-level 
and the p-level band. There our evidence for singlet superconductivity is very limited (at 
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least for frequencies larger than ~ 1/200). However, we have on that point found very clear 
evidence for superconductivity in the triplet sector. Following the procedure outlined above, 
we find consistently at small 'effective' temperature that any small terms in 0, in addition 
to the normal state solution, blow up at a rate which corresponds to a largest eigenvalue 
of ~ 1.8 of the matrix dF(iuj) n+l / dF(iuj) n . In fig. 3 we show the zero-temperature normal 
state solution for Ud = 4.5, /i = e p — €d = 4, calculated on a grid of points corresponding to 
an effective temperature of (3 = 200. The inset shows the most unstable eigenvector in the 



triplet sector. Its corresponding eigenvalue is \ max = 1-75 ||I3| . Superconducting order of 
this kind has been first proposed by Berezinskii |TJ|] in the context of 3 He, and, very recently 
by Coleman, Miranda and Tsvelik |D| for heavy-fermion superconductors. 

In conclusion, we have studied the infinite dimensional two-band Hubbard model with 
an exact diagonalization method, which has given very strong evidence for a superconduct- 
ing instability at low temperature. Given the excellent fit of the Green's functions, and 
the smallness of the finite size effect (dependence on n s ), it seems to us to be difficult to 
escape the conclusion that the two-band model is indeed superconducting, with a strongly 
frequency-dependent order parameter. We have taken every effort to check the programs 
(exact diagonalization and Lanczos) against each other and against the extrapolated Quan- 
tum Monte Carlo results at sufficiently high temperature. Close to the normal state solution 
we have also not been plagued by possible multiple solutions of the minimization in eq. (^j). 
In fact, we have also been able to find perfectly converged solutions of the fully super- 



conducting phase |n|, but there the extrapolation with n s , and the possible problem of 
multiple solutions are more critical. Therefore, the linear stability analysis presented here 
is the most convincing argument in favor of a superconducting instability of the two-band 
Hubbard model which we are able to give at the present time. 

Finally, as has been discussed elsewhere J?|, the method does for the time being not 
allow a precise calculation of the densities of states, and a renormalization group procedure is 
clearly called for. Such a method would allow the calculation of regular and superconducting 
densities of state. It would also allow calculations at finite temperature with much larger 
precision. The search for a numerical renormalization group procedure for D = oo is in our 
opinion an extremely exciting challenge to take up. 
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Figure Captions 

1. Zero temperature imaginary-time functions Gdo(iu) and G^ d {iu) (real parts: upper) 
(imaginary part: lower) vs u Ud = 8, // = e p — €d = 4 as calculated with the Lanczos 
algorithm with n s = 8 sites. Inset: Real part of the functions for n s = 4,5,6,7 and 
= 8. If G'dQ^uj) and Gdo(iu>) were exactly equal this would be the exact solution. 
Note the excellent fit of the two quantities for n s = 8, and the systematic improvement 
with increasing number of sites. A quasi-exact solution is obtained at n s =8. 

2. Largest eigenvector of the matrix dF(iu) n+1 /dF(iu) n close to the normal state solution 
of fig. 1. for n s = 6, 7 and 8 (singlet sector). The corresponding eigenvalues are 
A ma x ~ 2 in all three cases. 

3. Zero-temperature Green's function G{iu) in the normal state at Ud = 4.5, \i = e p —td = 
4 at n s = 7 (effective f3 = 200, n max = 64). The misfit between Gdo{iu>) and G°ftf(iu) 
is given by E = 1.5 x 10 -4 (!) (cf. eq. D) and the maximum difference between the 
two functions is 7 x 10~ 4 . The inset shows the most unstable eigenvector in the triplet 
sector at n s = 7. The corresponding eigenvalue is X ma x = 1-75. 
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